1_Dynamics/2_Time_domain_simulation/Numeric representation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The problem with simple math on computers\n",
    "\n",
    "Have you ever considered how computers store numbers? Can you explain why this happens?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999999"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a == 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = 0.125 + 0.125 + 0.125 + 0.125 + 0.125 + 0.125 + 0.125 + 0.125\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b == 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computers use base 2 instead of base 10\n",
    "\n",
    "You've heard that computers are all about ones and zeros, right? What does this mean?\n",
    "\n",
    "When I write a \"normal\" number like 123, what I mean is $1\\times10^2 + 2\\times 10^1 + 3 \\times 10^0$. This idea is called base 10 or decimal representation. Computers use binary or base 2 representation. This means you would write $101_2 = 5_10$, with the subscript representing the base. The math would work out as $1\\times 2^2 + 0\\times2^1 + 1\\times 2^0$, just like in the 123 example.\n",
    "\n",
    "This representation is exact for integers, but we run into problems when we use fractions. For instance, we all know that $1/3$ doesn't have a finite representation in decimals, since $1/3 = 0.\\overline{33} = 3\\times 10^{-1} + 3\\times10^{-2}+\\cdots$ forever. Notice that in base 3, 1/3 works out fine as $0.1_3$ since $1/3 = 1\\times3^{-1}$ exactly.  So here's the problem with writing $0.1$ in binary:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[This](http://bartaz.github.io/ieee754-visualization/) visualisation shows how IEEE floats are represented and indicates the repeating structure of the representation of 0.1.\n",
    "\n",
    "<img src=\"../../assets/0.1.png\">\n",
    "\n",
    "We can see that the binary representation is not finite, so the computer treats 1/10 more like a number like 1/7 (which we all know has an infinite decimal representation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a great deal more information on this issue at these pages:\n",
    "\n",
    "* [The Floating-Point Guide](http://floating-point-gui.de/) - this is an easy-to-read page with lots of examples\n",
    "* [What Every Computer Scientist Should Know About Floating-Point Arithmetic](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) - a more in-depth analysis of floating-point with lots of math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Built-in to Python\n",
    "\n",
    "The solution that Python supplies in the standard library is the decimal module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import decimal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = decimal.Decimal('0.1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 1/decimal.Decimal(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('1.0')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(a for i in range(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999999"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(0.1 for i in range(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sympy\n",
    "\n",
    "Sympy also has a solution in the form of a Rational object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = sympy.Rational('0.1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAC4AAAAUCAYAAADyWA/8AAAABHNCSVQICAgIfAhkiAAAAcFJREFUSInt1k2LjlEYB/DfmJK8PZjCgsI04xOIiIzFUFJGPoONbOQDjM1kJS8bK0rWhoWEYiGytHt4hJEixYIhLzVmFud6mO6ee7pv50nS/Ot03edc1/lf/859znUO/wkO4zzu4xOmcaXCvDWYwrkucK3DRbzBd0zgDFbONelxJJhEs0ayIxE7lMnVj3cRew2ncDf6T9BXNnEIA+jB7hrCb+I9ejO5bkXcscL46Ri/UEFL5WQN6ZdeyuTaFDEvsaDgW4bP+IIlOgT8CfZjIa5m8uwJexs/C75JPMBibKM7wkeklbiTybM5bKvE/yzsIPnCF2GftMe/ZXI1wn4s8bfHV5AvfBhLMZ7JUwU9YafJFz6CH7iRycPvFW2U+JfPjssR3osDUp0t+7118DTsYIl/IGyLPOG7pAuhW9vkXthhncvhDnzFIx0C6uCQVLauZ3DMxnOpFG7A0YLvpFS/L0sV7NeGb+NgNFiLvXghvTdIt+OJ+H4tvSN2lgipw9VGPx5itbQgTWyVbuEWtuNDp2Sj0qktaxMRtyX6x0tE1+EqYr10C7+VDv4rnMWqOXJVxlgk39gNsr+JpvT6m8c8/nXMAJtdflrva51YAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{1}{10}$$"
      ],
      "text/plain": [
       "1/10"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use `sympy.nsimplify`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = 1\n",
    "c = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = b/c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "float"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAC4AAAAUCAYAAADyWA/8AAAABHNCSVQICAgIfAhkiAAAAcFJREFUSInt1k2LjlEYB/DfmJK8PZjCgsI04xOIiIzFUFJGPoONbOQDjM1kJS8bK0rWhoWEYiGytHt4hJEixYIhLzVmFud6mO6ee7pv50nS/Ot03edc1/lf/859znUO/wkO4zzu4xOmcaXCvDWYwrkucK3DRbzBd0zgDFbONelxJJhEs0ayIxE7lMnVj3cRew2ncDf6T9BXNnEIA+jB7hrCb+I9ejO5bkXcscL46Ri/UEFL5WQN6ZdeyuTaFDEvsaDgW4bP+IIlOgT8CfZjIa5m8uwJexs/C75JPMBibKM7wkeklbiTybM5bKvE/yzsIPnCF2GftMe/ZXI1wn4s8bfHV5AvfBhLMZ7JUwU9YafJFz6CH7iRycPvFW2U+JfPjssR3osDUp0t+7118DTsYIl/IGyLPOG7pAuhW9vkXthhncvhDnzFIx0C6uCQVLauZ3DMxnOpFG7A0YLvpFS/L0sV7NeGb+NgNFiLvXghvTdIt+OJ+H4tvSN2lgipw9VGPx5itbQgTWyVbuEWtuNDp2Sj0qktaxMRtyX6x0tE1+EqYr10C7+VDv4rnMWqOXJVxlgk39gNsr+JpvT6m8c8/nXMAJtdflrva51YAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{1}{10}$$"
      ],
      "text/plain": [
       "1/10"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.nsimplify(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why isn't math always done in base 10?\n",
    "\n",
    "The extra precision comes at a cost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5 ms ± 321 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "a = 0.1\n",
    "s = 0\n",
    "for i in range(100000):\n",
    "    s += a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.8 ms ± 246 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "a = decimal.Decimal('0.1')\n",
    "s = decimal.Decimal(0)\n",
    "for i in range(100000):\n",
    "    s += a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.03 s ± 166 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "a = sympy.Rational(1, 10)\n",
    "s = 0\n",
    "for i in range(100000):\n",
    "    s += a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using sympy rationals is about a thousand times slower than using built-in Python `float`s."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forcing rounding of exact representations\n",
    "\n",
    "If an equation results in an \"Exact\" answer which isn't \"useful\", like $\\sqrt{3}x$, we can approximate that using `sympy.N`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sympy.Symbol('x')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACsAAAAbCAYAAADlJ3ZtAAAABHNCSVQICAgIfAhkiAAAAiZJREFUWIXt1k2ITWEYB/DfzCBhFgwNNWXB1JCMZCGNFfKxw6SUha2yEFZSZkGXpklNPmKh2VCslAWpWWmklCEiiqamZhAGCx/jc/G+073d7j1zjrm3sZh/ne57nuf59/+f57zneS9TqA5qMtbX4/Mk6mdCOzZXU6CSOIMZkyVem6G2Ll6jVfIyLrKYXYe71TJSaeTQMJkGsnS2Ae+rZSQN0pptxvMyuVPoxSC+4gP6cUyF30TaOXcIN/CiRG4UD/AUbzEba7EGQ3E9OGGnGXAhITezTPwE/uBcpUzUFq2PY2tRzTzJe/Vbmfi1+NucwL0tPNCOongNemLuZCniFmHvXSqK70FbgmA5HI1iXQk1rfglbKG6gnhX5F5MElgl7LvCjp+V7kM8jA6cxp0o9ggLxuH1xNq98f5IvL+aRndAvpPT0Z3CKLyOImPXTTSm4DUJU2QA+yP3lpTHejc643oTdqY0O4ZGbBdG3RBWp+Dk5B+yD7PSim2QH1GdmJPFaQEW4zuepKg9KG+2JYvINIxgmYmPnf5oYH5CzW78xnCsPZ9V5DKuYN8/GCzEm2hgbpn8NuFQeSx8iM/wQ8bu7ooiTePUtWBhiXit/KHQV4bbhi94hUUx1h4517OYrZfu7+ABoRO9wkzMCXP6ZRQdxvISvFZ8jPklRbn7kbs+i+Fyr64QK4Q5/BDv8BOfomCHcPoVY6kw5kawskR+YzR7L4vZKUzhf8RfH812nZJaPpEAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\sqrt{3} x$$"
      ],
      "text/plain": [
       "√3⋅x"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr = sympy.sqrt(3)*x\n",
    "expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACMAAAAWCAYAAABKbiVHAAAABHNCSVQICAgIfAhkiAAAAcRJREFUSInt1TtoFFEUxvFfIohoLHwRBcHCIEGiAbGw0EasBAvFxk57CxEbRXAbiSJBEIxgIWlTKVj4IpUEhIBRFMVCEQImUYmPQvGxrsU9yw7jbMi6G1OYDy7MnHPmu/95nDP8xzqJUXzGO9xEz3zB3MGRANiC65jEyvkCyqoDZeyD9lzyPIYxjq+YxhjOYNUcwCwPhmloyyW/4yGe4S2WYQe2400cj7cQZgibwr+cTy6pc9FZVDDQQpALmEBXoxf2Bsy9XPxuxA/k4m0YjNy5Ar9+TGFzoyBwOoz7CyDL0mtdlNusgqsFXpcaBTmBEi7ifhg/xpqC2sHIH47zU3E+5M8GGZBmzG6szayOmWAmw7C6bqGzTu16qfNe42jU38bigtpKnVWaCaaqTuzHC6mbttWp68sYj2DpbMz/VhvwDU/r5I9nYLrnEqSqsdhsdS5+CL+kNq3gyr+AmYrNVmRie6VB+UT6uJ/jhxY8nW7py86rXW3ojWTiO/EFr7AuYgej7kazMMekuxqW5kMfruFlbDChNht68TFiG3M+o1G/qxmYHlzGI7zHT3wK85Lab75Lav0P2FrgsydgHjQDs6AFzUa/AZZEdZ+WkX7QAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$3 x^{2}$$"
      ],
      "text/plain": [
       "   2\n",
       "3⋅x "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.simplify(expr**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADcAAAAWCAYAAABkKwTVAAAABHNCSVQICAgIfAhkiAAAArxJREFUWIXt10+IVXUUB/CPRhL4j0xzAsPwqQz5ZyAkhGphhQuhjbhxV+hGKBRqowjOJlSm6Y+WmIuSdq4yWqSFbiIRBjQR1BGsCaVxMi2llKZkXPzOZZ733eu775mMi/nCj/vuOd9zfufc+zvn3Mc4HkpsRh9u4Aq+xuIxjeh/xGG8ISW0BF/iMmaMZVAPClNwG69lgok5wk4cwUXcwjWcxDY80caGc/AZfsU/GMCHeLwNX80wVcrnWiaYkCMM4wTO4DdMxnIsiwCXS4lXQQ3H8CS+wjk8jxXoxwu42l4ehTiAhRHr7SLCYyWG72IEe1rY7HDYvJWTvx/yvS34aoYeDGJ+O8ZdEdB3Ffnzgv+zxqM/FX/hb+lkZPg2bFbn+BOwP3Q7CvbqxRCerRhbA7aG896K/PXB/7REn73VV+pkXdJxOoNH6uS9wd1X4GeXNhJ7B934AN+H81OYVdG+J2zeLtF/HPoNOfn+kL8e91vi/oDGE7BHmnEvo6NuTWkW3OVwmq1vMLuZUR32hd36En1Ww5tz8jlSlx7Am8E5hEkFPkZKVndGyD+NDB3SWe+QamCeNBKeu1dGLSDr0iM5+SVpVMzFbqnbrpa6eJGPotWdEcqSyzAkTf6V0pz7omLw1+M6vUQ/Lcerx5W63+tws+KeDWiWXIZfpEJfhJkV+P1xXViiXxDX8zn5WrwnlQVsrBjffWNIOkZVvi5qmo+Cm+4eBauk43daalxn8S867yvqQKdUY3lMNNoAfsjpamH3aIFdK0P8RSnZn/BUyNYE72DlDO6BTdKTOiJ1u+3Sd+GF2GRQ4zwZCN0zBf5qRt/2wfB3NO77jX6rduHP8F/L+egL/kvtJpVhMT7Bj/gd/0kF3yd1oKK/EgPKk4On8bkU+LBUux/V+Zov1dcfWFpg/2r4P95aKuMYxzjGAncAsEKuIt6sOvgAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$3.0 x^{2}$$"
      ],
      "text/plain": [
       "     2\n",
       "3.0⋅x "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.simplify(sympy.N(expr, 3)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}